#!/usr/bin/env python
"""
Runs Monte Carlo program for the case of a triangular lattice quantum
transverse field ising model with L**2 spins by mapping it to the
classical 3d ising model with Z layers. The interaction energy and the
field energy are measured seperately and stored.
"""


from __future__ import division
import numpy as np
import numpy.random as ra
import math
from datetime import datetime

def compute_lattice_energy(S, L, Z):
    EB = 0
    EJ = 0
    for k in range(Z):
        for i in range(L):
            for j in range(L):
                EJ += - S[k, i, j] * sum_h_neigh(S, k, i, j)
                EB += - math.exp(- 2 * lamb * S[k, i, j] * S[u[k], i, j])
    EJ *= 0.5 * J_orig
    EB *= B 
    return EJ, EB
    

def sum_h_neigh(S, k, i, j):
    h = S[k, n[i], j] + S[k, s[i], j] + S[k, i, e[j]] + S[k, i, w[j]] + S[k, n[i], w[j]] + S[k, s[i], e[j]]
    return h

def sum_v_neigh(S, k, i, j):
    v = S[u[k], i, j] + S[d[k], i, j]
    return v


t_start = str(datetime.now()) # Stores the time when program starts. Used for measuring running time.
L = 4 # Total no. of spins is L^2
tot_sweeps = 10000
J_orig = 0 # Decides whether the system is ferromagnetic (+1) or antiferromagnetic (-1)
B = 0.2 # Magnetic field
dt = 1/8 # Decides magnitude of error

# Periodic boundary
n = [L-1] + range(L-1)
w = n
s = range(1, L) + [0]
e = s



temp_range = np.arange(0.5, 0.6, 0.2)
prefix = 'time_series_f_tri_L_'+`L`+'_dt_%.3f_ising_T_'%(dt)
t_trans = 1000

for T in temp_range:
    Z = int(1/(T * dt)) # Determine the number of layers
    print Z
    # Define neigbhours for the vertical layers
    u = [Z - 1]  + range(Z - 1)
    d = range(1, Z) + [0]
    
    S = np.sign((ra.uniform(0, 1, [Z, L, L]) - 0.5)) # Initialize the lattice
    EJ = 0 #Energy in spins
    EB = 0 #Energy in field
    accepted = 0 # No. of accepted changes

    lamb = - 0.5 * math.log( math.tanh( B * dt ) ) # define \lambda

    EJ, EB = compute_lattice_energy(S, L, Z)
    print EJ, EB
    EJ_temp = []
    EB_temp = []
    M_temp = []
    sweep = 0
    while (sweep < tot_sweeps):
        for k in range(Z):
            for i in range(L):
                for j in range(L):
                    dEJ = 2 * S[k, i, j] * (J_orig * sum_h_neigh(S,k,i,j))
                    dS = 2*S[k, i, j]*( lamb * sum_v_neigh(S, k, i, j) )  + dEJ * dt
                    if ( ra.uniform(0, 1) < math.exp(-dS)/(1.0 + math.exp(-dS) ) ):
                        accepted += 1
                        dEB = - 2 * B * ( math.sinh(2*lamb*S[k, i, j]*S[u[k], i, j]) + math.sinh(2*lamb*S[k, i, j]*S[d[k], i, j]) )
                        S[k, i, j] *= -1
                        EJ += dEJ
                        EB += dEB
        #EJ1, EB1 = compute_lattice_energy(S, L, Z)
        if (sweep > t_trans):
            M_temp.append(np.sum(S))
            EJ_temp.append(EJ)
            EB_temp.append(EB)
        sweep += 1

    EJ_temp = np.array(EJ_temp)
    fname = '../3d_tri/data/EJ_' + prefix + '%.2f_B_%.2f.data'%(T, B)
    EJ_temp = EJ_temp/(L*L*Z) # Calculate energy per site
    EJ_temp.tofile(fname, sep=', ')
    EB_temp = np.array(EB_temp)
    EB_temp = EB_temp/(L*L*Z) # Calculate energy per site
    fname = '../3d_tri/data/EB_' + prefix + '%.2f_B_%.2f.data'%(T, B)
    EB_temp.tofile(fname, sep=', ')
    M_temp = np.array(M_temp)
    M_temp = M_temp/(L*L*Z) # Calculate energy per site
    fname = '../3d_tri/data/M_' + prefix + '%.2f_B_%.2f.data'%(T, B)
    M_temp.tofile(fname, sep=', ')

print np.mean(EB_temp)
t_end = str(datetime.now())
fname = '/home/vikram/' + prefix +'.t'
f = open(fname,'w')
info = t_start + '\n' + t_end
f.write(info)
f.close()
